# GRAPHS 1-6 FOR CIRCULAR DATA ARTICLE

setwd("C:/H/circular data/graphs")
library(arm); library(circular); library(MASS)

# graph1.bw.ps
pitt <- read.dta("C:/H/circular data/guncrime/02895-0005-Data.dta")
pitt.hour <- circular(2*pi*pitt$HOUR/23,type="angles",units="radians",rotation="clock")
postscript("graph1.bw.ps")
par(mar=c(2,2,1,1),mfrow=c(1,2))
bar.out <- barplot(table(pitt$HOUR), space=0,col="lightgray",xlab="",ylab="")
set.seed(668)
plot(sample(pitt.hour,800), stack=FALSE, shrink=1.3, cex=1.03,
    axes=FALSE,tol=0.8,zero=c(rad(90)),bins=24,ylim=c(0,3))
points(sample(pitt.hour,800), rotation = "clock", zero =c(rad(90)), 
    col = "1", cex=1.03, stack=TRUE )    
text(0,0.8,"24"); text(0,-0.8,"12"); text(0.8,0,"6"); text(-0.8,0,"18")   
dev.off()

# graph2.bw.ps
pitt.hour.short <- as.circular(sample(pitt.hour,size=1500))
postscript("graph2.bw.ps",width=5.2,height=5.2)
par(mfrow=c(1,1),mar=c(4,4,4,4))
plot(pitt.hour.short-pi/2+0.005+rnorm(length(pitt.hour.short),0,0.02),pch=3,col="darkgrey")
mtext(side=3,line=0.65,"Midnight",cex=0.9); mtext(side=1,line=0.65,"Noon",cex=0.9)
mtext(side=2,line=0.65,"9PM",cex=0.9); mtext(side=4,line=0.65,"3AM",cex=0.9);
par(new=TRUE,col.axis="white")
rose.diag(pitt.hour-pi/2+rnorm(length(pitt.hour),0,0.02),bins=24,shrink=0.33,xlim=c(-2,2),ylim=c(-2,2),axes=FALSE,prop=1.5)
dev.off()




# graph3.bw.ps
postscript("graph3.bw.ps")
circ.dens <- density(pitt.hour+3*pi/2, bw=10)
par(mar=c(2,4,0,1),mfrow=c(1,1))
# plot(pitt.hour, stack=FALSE, shrink=1.5, cex=1.03,axes=FALSE,tol=0.8,zero=c(rad(90)),bins=24,ylim=c(-1,1))
plot(sample(pitt.hour,1), stack=TRUE, shrink=.35, cex=0, sep=0.0,
    axes=FALSE,tol=.8,zero=c(0),bins=24,ylim=c(-1.29,3), ticks=TRUE, tcl=.075)
lines(circ.dens, col="darkgrey", lwd=3)
text(0,0.8,"24", cex=2); text(0,-0.8,"12",cex=2); text(0.8,0,"6",cex=2); text(-0.8,0,"18",cex=2)
dev.off() 



# graph4.bw.ps
postscript("graph4.bw.ps")
par(mar=c(1,1,1,1),mfrow=c(1,3),cex.lab=2)
# 1st plot
circ.draw1 <- rwrappednormal(300,mu=circular(0),sd=1)
circ.dens1 <- density(circ.draw1, bw=5)
plot(circ.dens1, xlim=c(-1.5,1.5), template = "clock", 
    rotation = "clock", col="darkgrey", lwd=3, zero =c(rad(90)), main=" ")
points(circ.draw1, rotation = "clock", zero =c(rad(90)), col = "1", pch=6, cex=1.03 )
mtext(side=1,line=-31.7,expression(mu==0),cex=1.2)
mtext(side=1,line=-27.7,expression(sigma^2==1),cex=1.2)
# 2nd plot
circ.draw1 <- rwrappednormal(300,mu=circular(pi/4),sd=1/3)
circ.dens1 <- density(circ.draw1, bw=5)
plot(circ.dens1, xlim=c(-1.5,1.5), template = "clock", 
    rotation = "clock", col="darkgrey", lwd=3, zero =c(rad(90)), main=" ")
points(circ.draw1, rotation = "clock", zero =c(rad(90)), col = "1", pch=6, cex=1.03 )
mtext(side=1,line=-31.7,expression(mu==pi/4),cex=1.2)
mtext(side=1,line=-27.7,expression(sigma^2==1/3),cex=1.2)
# 3rd plot
circ.draw1 <- rwrappednormal(300,mu=circular(5*pi/4),sd=0.75)
circ.dens1 <- density(circ.draw1, bw=5)
plot(circ.dens1, xlim=c(-1.5,1.5), template = "clock", 
    rotation = "clock", col="darkgrey", lwd=3, zero =c(rad(90)), main=" ")
points(circ.draw1, rotation = "clock", zero =c(rad(90)), col = "1", pch=6, cex=1.03 )
mtext(side=1,line=-31.7,expression(mu==5*pi/4),cex=1.2)
mtext(side=1,line=-27.7,expression(sigma^2==3/4),cex=1.2)
dev.off()

# graph5.bw.ps
postscript("graph5.bw.ps")
set.seed(666)
circ.draw100 <- rwrappedcauchy(100, mu=circular(0))
circ.draw10k <- rwrappedcauchy(10000, mu=circular(0))
circ.dens <- density(circ.draw10k, bw=30)
par(mar=c(1,1,1,1),mfrow=c(1,1),cex.lab=2)
plot(circ.dens, xlim=c(-1.5,1.5), template = "clock", shrink=1.2,
    rotation = "clock", col="darkgray", lwd=3, zero =c(rad(90)), main=" ")
points(circ.draw100, rotation = "clock", zero =c(rad(90)), col = "1", pch=1, cex=1 )
dev.off()


# graph6.bw.ps
postscript("C:/H/circular data/graphs/graph6.bw.ps", paper="letter")
par(mar=2*c(.9,2.2,.1,1),mfrow=c(3,3),cex.lab=2)
N <- 1000
X <- (runif(N)-.5) * 2
X.d <-c(X,X) # duplicate X
# generate data 1
Y1 <- as.circular(matrix(NA,N,1), type=c("angles"), units=c("radians"),
    rotation = c("clock"), template="geographics")
for(i in 1:N) {
    Y1[i] <- rvonmises(1, mu=circular(0 + 2*atan(0 * X[i])), kappa = exp(-1) ,
                 control.circular = list(units = "radians"))
 }
# generate data 2 
Y2 <- as.circular(matrix(NA,N,1), type=c("angles"), units=c("radians"),
    rotation = c("clock"), template="geographics")
for(i in 1:N) {
    Y2[i] <- rvonmises(1, mu=circular(1 + 2*atan(2 * X[i])), kappa = 2 ,
                 control.circular = list(units = "radians"))
 } 
# generate data 3 
Y3 <- as.circular(matrix(NA,N,1), type=c("angles"), units=c("radians"),
    rotation = c("clock"), template="geographics")
for(i in 1:N) {
    Y3[i] <- rvonmises(1, mu=circular(0 + 2*atan(4 * X[i])), kappa = 2 ,
                 control.circular = list(units = "radians"))
 } 
#scatterplot 1
Y1.d  <- c(Y1,Y1 + 2*pi) 
plot(X.d, deg(Y1.d), ylim=c(0,730), xlab="x", ylab="Degrees", pch=".") 
#scatterplot 1
Y2.d  <- c(Y2,Y2 + 2*pi) 
plot(X.d, deg(Y2.d), ylim=c(0,730), xlab="x", ylab=" ", pch=".") 
#scatterplot 1
Y3.d  <- c(Y3,Y3 + 2*pi) 
plot(X.d, deg(Y3.d), ylim=c(0,730), xlab="x", ylab=" ", pch=".") 
# density 1
res25 <- density(Y1, bw=25) 
plot(res25, points.plot=FALSE, c(-1.15,1.15), template = "clock", shrink=1.22,
    rotation = "clock", col="darkgrey", lwd=3, zero =c(rad(90)), main=" ", ylab=" ") 
mtext(side=3,line=-6,expression(mu==0),cex=2)
mtext(side=3,line=-10,expression(beta==0),cex=2)
mtext(side=3,line=-14,expression(kappa==e^-1),cex=2)
# density 2   
res25 <- density(Y2, bw=25) 
plot(res25, points.plot=FALSE, c(-1.15,1.15), template = "clock", shrink=1.22,
    rotation = "clock", col="darkgrey", lwd=3, zero =c(rad(90)), main = " ", ylab=" ") 
mtext(side=3,line=-6,expression(mu==1),cex=2)
mtext(side=3,line=-10,expression(beta==2),cex=2)
mtext(side=3,line=-14,expression(kappa==2),cex=2)    
# density 3          
res25 <- density(Y3, bw=25) 
plot(res25, points.plot=FALSE, xlim=c(-1.15,1.15), template = "clock", shrink=1.22,
    rotation = "clock", col="darkgrey", lwd=3, zero =c(rad(90)), main = " ", ylab=" ") 
mtext(side=3,line=-6,expression(mu==0),cex=2)
mtext(side=3,line=-10,expression(beta==4),cex=2)
mtext(side=3,line=-14,expression(kappa==2),cex=2) 
# R-loglik 1    
beta.trial <- seq(-7,7,length.out=1000)
R <- matrix(NA,length(beta.trial),1)
for(i in 1:length(beta.trial))   { 
 S <- sum(sin(Y1 - 2 * atan(X * beta.trial[i])))/N
 C <- sum(cos(Y1 - 2 * atan(X * beta.trial[i])))/N
 R[i] <- sqrt(S^2 + C^2)
 R
}
plot(beta.trial,R, type="l", xlab= " ", ylab="R Component of LL") 
M1 <- circ.lin.reg(X,Y1,c(0), trace=T, tol = 1e-5)
abline(v=M1$beta, col="darkgrey", lwd=3, lty=2) 
# mtext(side=1,line=+1,expression(beta),cex=.9)   
# R-loglik 2    
R <- matrix(NA,length(beta.trial),1)
for(i in 1:length(beta.trial))   { 
 S <- sum(sin(Y2 - 2 * atan(X * beta.trial[i])))/N
 C <- sum(cos(Y2 - 2 * atan(X * beta.trial[i])))/N
 R[i] <- sqrt(S^2 + C^2)
 R
}
plot(beta.trial,R, type="l", xlab= " " , ylab=" ") 
M2 <- circ.lin.reg(X,Y2,c(2), trace=T, tol = 1e-5)
abline(v=M2$beta, col="darkgrey", lwd=3, lty=2)  
# mtext(side=1,line=+1,expression(beta),cex=.9)     
# R-loglik 3     
R <- matrix(NA,length(beta.trial),1)
for(i in 1:length(beta.trial))   { 
 S <- sum(sin(Y3 - 2 * atan(X * beta.trial[i])))/N
 C <- sum(cos(Y3 - 2 * atan(X * beta.trial[i])))/N
 R[i] <- sqrt(S^2 + C^2)
 R
}
plot(beta.trial,R, type="l", xlab= " " , ylab=" ")  
M3 <- circ.lin.reg(X,Y3,c(4), trace=T, tol = 1e-5)   
abline(v=M3$beta, col="darkgrey", lwd=3, lty=2) 
# mtext(side=1,line=+1,expression(beta),cex=.9)   
    
dev.off()
# circ.lin.reg

library(circular)


circ.lin.reg <- function(x, theta, beta0, trace = FALSE, print = TRUE, tol = 1e-10, maxiter=1000)
{
# THIS CODE IS HEAVILY BASED ON S-CODE FROM ULRIC LUND'S CIRCSTATS V2.0 (2004)
        if(is.vector(x))
                x <- cbind(x)
        n <- length(theta)
        betaPrev <- coef(lm(theta~x))[2:(dim(x)[2]+1)]
        # CHANGE 1:1 TO 1:4 FOR DIFFERENT STARTING VALUES
    for(i in 1:1)  {
    if(i==2) betaPrev <- betaPrev
    if(i==3) betaPrev <- betaPrev+1
    if(i==4) betaPrev <- betaPrev-1   
    if(i==1) betaPrev <- rep(0,length(betaPrev)) 
        S <- sum(sin(theta - 2 * atan(x %*% betaPrev)))/n
        C <- sum(cos(theta - 2 * atan(x %*% betaPrev)))/n
        R <- sqrt(S^2 + C^2)
         mu <- atan2(S, C)
        k <- A1inv(R)
        diff <- tol + 1
        iter <- 0
        S.function <- function(betaPrev,x) {
                2/(1 + (t(betaPrev) %*% x)^2) }
        while(diff > tol & iter < maxiter) {
                iter <- iter + 1
                u <- k * sin(theta - mu - 2 * atan(x %*% betaPrev))
                A <- diag(k * A1(k), nrow = n)

                g.p <- diag(apply(x, 1, S.function, betaPrev = betaPrev), nrow = n)
                D <- g.p %*% x
                betaNew <- lm(t(D) %*% (u + A %*% D %*% betaPrev) ~ t(D) %*% A %*% D - 1)$coefficients
                diff <- abs(max(betaNew - betaPrev))
                breaked = 0
                if(iter==1000) {
                breaked = 1
                }
                if(is.na(diff)==TRUE) {
                breaked = 1
                break
                }
                if(max(betaNew) > 100)  {
                breaked = 1
                break
                }
                betaPrev <- betaNew
                S <- sum(sin(theta - 2 * atan(x %*% betaPrev)))/n
                C <- sum(cos(theta - 2 * atan(x %*% betaPrev)))/n
                R <- sqrt(S^2 + C^2)
                 mu <- atan2(S, C)
                # mu <- asin(S/R)
                k <- A1inv(R)
                if(trace == T) {
                        log.lik <-  - n * log(2*pi*I.0(k)) + k * sum(cos(theta -
                                mu - 2 * atan(x %*% betaNew)))
                        cat("Iteration ", iter, ":    Log-Likelihood = ",
                                log.lik,"mu ", mu, "k ", k, "b ", betaNew, "\n")
                        cat("Starting values ", i, "\n")
                }
        }
        log.lik <-  - n * log(2*pi*I.0(k)) + k * sum(cos(theta - mu - 2 * atan(
                x %*% betaNew)))
        log.lik.old <-  - n * log(2*pi*I.0(k)) + k * sum(cos(theta - mu - 2 * atan(
                x %*% betaPrev)))
        cov.beta <- solve(t(D) %*% A %*% D)
        se.beta <- sqrt(diag(cov.beta))
        se.kappa <- sqrt(1/(n * (1 - A1(k)^2 - A1(k)/k)))
        circ.se.mu <- 1/sqrt((n - ncol(x)) * k * A1(k))
        z.values <- abs(betaNew/se.beta)
        p.values <- (1 - pnorm(z.values))*2
        result.matrix <- cbind(Coef = betaNew, SE = se.beta, Z = z.values,
                p = p.values)
        dimnames(result.matrix) <- list(dimnames(x)[[2]], c("Coef", "SE", "|z|",
                "p"))
        cat("\n", "Circular-Linear Regression", "\n", "\n")
        print(result.matrix)
        cat("\n", "\n")
        betaNew <- as.matrix(betaNew)
        dimnames(betaNew) <- list(dimnames(x)[[2]], c("Estimate"))
        list(mu = mu, kappa = k, beta = betaNew, log.lik = log.lik, log.lik.old = log.lik.old, 
                circ.se.mu = circ.se.mu, se.kappa = se.kappa, cov.beta = cov.beta,
                se.beta = se.beta, result.matrix = result.matrix, breaked=breaked)
   }
}
####
rm(list=ls())
options(digit=4)
library(circular)
library(MASS)
library(arm)
library(mice)
library(Design)

#library(Cairo)


violence <- read.dta("C:/H/circular data/violence/violence2.dta")
summary(violence)
# convert angles
violence$direction.rad <- violence$YEARDAY / 365 * 2* pi
violence$direction.rad <- as.circular(violence$direction.rad, units = "radians", 
    rotation = "clock", zero =c(rad(0)), type="angles")
summary(violence$direction.rad)

source("C:/H/circular data/circ.lin.reg.R")

viol.short <- violence[ ,c("YEARDAY","direction.rad","NATURE_OF_THE_TARGET","NUMBER_OF_ATTACKERS",
                    "MOTIVATION_FOR_ATTACK","NUMBER_KILLED_TARGET","NUMBER_KILLED_ATTACKER",
                    "TYPE_GROUP_TARGET","TYPE_GROUP_ATTACKER","YEAR_OF_ISSUE","GROUP_VIOLENCE",
                    "PROPERTY_AS_TARGET","NEWSPAPER_IDENTIFICATION","CODER_IDENTIFICATION")]
                    
# do not run again
#viol.imp <- mice(viol.short, m = 100, seed=666)


sink("C:/H/circular data/violence/result.txt")

imps <- 50
coefs <- matrix(NA,12,imps)
se.beta <- matrix(NA,12,imps)
log.lik <- matrix(NA,1,imps)
log.lik.old <- matrix(NA,1,imps)
mu <- matrix(NA,1,imps)
kappa <- matrix(NA,1,imps)
circ.se.mu <- matrix(NA,1,imps)
se.kappa <- matrix(NA,1,imps)
breaked <- matrix(NA,1,imps)

for (i in 1:imps)    {
viol.data <- complete(viol.imp, action=i)

year <- viol.data$YEAR_OF_ISSUE


N <- length(year)
Y <- viol.data$direction.rad

year2 <- year-min(year)

killed.target <- viol.data$NUMBER_KILLED_TARGET
killed.target[killed.target==97] <- 0
killed.target[killed.target==51] <- 75
killed.target[killed.target==52] <- 300
killed.attack <- viol.data$NUMBER_KILLED_ATTACKER
killed.attack[killed.attack==97] <- 0
killed.attack[killed.attack==51] <- 75
killed.attack[killed.attack==52] <- 300

killed.target.bin <- ifelse(killed.target >= 1,1,0)
killed.attack.bin <- ifelse(killed.attack >= 1,1,0)

property.target <- viol.data$PROPERTY_AS_TARGET
property.non <- ifelse(property.target == 0,1,0)
property.per <- ifelse(property.target == 1,1,0)
property.bus <- ifelse(property.target == 2,1,0)
property.gov <- ifelse(property.target == 3,1,0)
property.rel <- ifelse(property.target == 4,1,0)
property.oth <- ifelse(property.target == 7,1,0)


attack.type <- viol.data$TYPE_GROUP_ATTACKER

attack.rac <- ifelse(attack.type == 0,1,0)
attack.lab <- ifelse(attack.type == 1,1,0)
attack.oth <- ifelse(attack.type == 2 | attack.type == 3 | attack.type == 6 | attack.type == 7,1,0)
attack.pol <- ifelse(attack.type == 4,1,0)
attack.soc <- ifelse(attack.type == 5,1,0)
#attack.off <- ifelse(attack.type == 6,1,0)
#attack.oth <- ifelse(attack.type == 7,1,0)

                               
target.type <- viol.data$TYPE_GROUP_TARGET

target.rac <- ifelse(target.type == 0,1,0)
target.lab <- ifelse(target.type == 1,1,0)
target.bus <- ifelse(target.type == 2,1,0)
target.oth <- ifelse(target.type == 3 | target.type == 5 | target.type == 7,1,0)
target.pol <- ifelse(target.type == 4,1,0)
# target.soc <- ifelse(target.type == 5,1,0)
target.off <- ifelse(target.type == 6,1,0)
#attack.oth <- ifelse(attack.type == 7,1,0)   


number.attack <- viol.data$NUMBER_OF_ATTACKERS
num.attack.bin <- ifelse(number.attack >= 2,1,0)
num.attack.bin10 <- ifelse(number.attack >= 7,1,0)


###
X <- cbind(year2/100, killed.target.bin, num.attack.bin10,
    attack.rac,attack.lab,           attack.oth,attack.soc,
    target.rac,target.lab,target.bus,target.oth,target.off)  


# baseline: attack.pol, target.pol

out1<-circ.lin.reg(X,Y, trace=T, tol = 1e-5, maxiter=100)
log.lik[i] <- out1$log.lik
log.lik.old[i] <- out1$log.lik.old
coefs[,i] <- out1$beta
se.beta[,i] <- out1$se.beta
mu[i] <- out1$mu
kappa[i] <- out1$k
circ.se.mu[i] <- out1$circ.se.mu
se.kappa[i] <- out1$se.kappa
breaked[i] <- out1$breaked
}
sink()

index <- c(4,6,8,9,11,15,25,29,30,31,33,37,38,40,43,47,49,50) 

index.mu <- index[(mu[index] < (-2))]
index.ll <- index.mu[(log.lik[index.mu] > (-1337))]
index.ll
index.new <- index.ll
coefs[,index.new]
mu[index.new]
log.lik[index.new]

out.null<-circ.lin.reg(matrix(1,N,1),Y, trace=T, tol = 1e-10, maxiter=100)


mean(log.lik[index.new]) # MEAN LL FULL MODEL
out.null$log.lik # LL NULL MODEL
(dev1 = 2 * (mean(log.lik[index.new]) - out.null$log.lik)) # Fahrmeier/Tutz DEVIANCE
(dev2 = 2 * (-mean(log.lik[index.new]))) # RESIDUAL DEVIANCE
(dev3 = 2 * (0- out.null$log.lik)) # NULL DEVIANCE
(AIC <- 2*dim(X)[2] - 2*mean(log.lik[index.new])) # AIC

pars <- rbind(mu[index.new],coefs[,index.new],kappa[index.new])
pars.se <- rbind(circ.se.mu[index.new],se.beta[,index.new],se.kappa[index.new])

# USE THIS FUNCTION TO GET THREE REQUIRED VECTORS:
impute.coef.vec <- apply(pars,1,mean)

between.var <- apply(pars,1,var)

within.var <- apply(pars.se^2,1,mean)


# COMPUTE THE STANDARD ERROR VECTOR
m <- length(index.new)
impute.se.vec <- sqrt(within.var + ((m+1)/m)*between.var)

# THE DEGREES OF FREEDOM FOR THE T-STATISTIC NEEDS TO BE ADJUSTED.
# SEE LITTLE AND RUBIN (1987), PAGE 257
impute.df <- (m-1)*(1 + (1/(m+1)) * within.var/between.var)^2 

# REGRESSION TABLE:
out.table <- round( cbind( impute.coef.vec,impute.se.vec,impute.coef.vec/impute.se.vec,
        (1-pt(abs(impute.coef.vec/impute.se.vec),impute.df))*2), 3 )
out.table
dimnames(out.table) <- list( c("(mu)","year","some killed","large group",
        "attack: race","attack: labor","attack: other","attack: social",
        "target: race","target: labor","target: business","target: other","target: officials",
        "kappa"),
        c("Estimate","Std. Error","t value","Pr(>|t|)") ) 
        
        
# CODE TO GENERATE FIGURE

# MI
viol.imp <- mice(viol.short2, m = 5, seed=666)
viol.data <- complete(viol.imp, action=1)

# BIG TABLE JUST FOR EVALUATION
viol.imp <- read.dta("./Article.Circular/ICPSR.10811281/viol.data.imp.dta")
categories <- 8
( count.table <- table(viol.imp$TYPE_GROUP_TARGET, viol.imp$TYPE_GROUP_ATTACKER)[categories:1,] )

# MAKE A NICE GRAPH FOR CATEGORICAL DATA
attack.type <- viol.imp$TYPE_GROUP_ATTACKER
attack.type[attack.type %in% c(3,6,7)] <- 2
target.type <- viol.imp$TYPE_GROUP_TARGET
target.type[target.type %in% c(5,7)] <- 3
( count.table <- table(target.type, attack.type)[6:1,] )
cat.rows <- nrow(count.table); cat.cols <- ncol(count.table)

postscript("./Article.Circular/terror.type.ps")
par(mar=c(3,3,2,3),cex=1.2)
plot(c(0,cat.cols),c(0,cat.rows),type="n",xlim=c(0,cat.cols+0.5),ylim=c(0,cat.rows+0.5),
    bty="n",xaxt="n",yaxt="n",xlab="",ylab="")
mtext(side=1,cex=1.3,line=1,"Attack Group Type"); mtext(side=2,cex=1.3,line=1,"Target Group Type")
text(0.5,6.25,"RAC"); text(5.25,0.5,"RAC")
text(1.5,6.25,"LAB"); text(5.25,1.5,"LAB")
text(2.5,6.25,"OTH"); text(5.25,2.5,"BUS")
text(3.5,6.25,"POL"); text(5.25,3.5,"OTH")
text(4.5,6.25,"SOC"); text(5.25,4.5,"POL")
                      text(5.25,5.5,"OFF")
num.lines <- 500
for (i in 1:cat.cols)  {
    for (j in 1:cat.rows)  {
        for (k in 1:num.lines) segments( (i-1)+k/num.lines, (j-1), (i-1)+k/num.lines, j, 
            col=colors()[357-count.table[(cat.rows-j+1),i]])
    text( (i-1)+.5, (j-1)+.5, count.table[(cat.rows-j+1),i], cex=1.2 )
    }
}
dev.off()        
        
####
rm(list=ls())
options(digit=4)
library(circular)
library(MASS)
library(arm)
#library(Cairo)
germany <- read.dta("C:/H/circular data/directional/germandata.dta")
attach(germany)  
summary(germany)
# convert counterclockwise angles
germany$direction.rad <- (360-germany$direction) / 360 * 2* pi
germany$direction.rad <- as.circular(germany$direction.rad, units = "radians", 
    rotation = "clock", zero =c(rad(0)), type="angles")
summary(germany$direction.rad)
detach(germany)
attach(germany)
germany$direction.rad 


# BLACK AND WHITE
#### graph density german bundestag
set.seed(668)
res5 <- density(direction.rad[pds==1], bw=5)
plot(res5, points.plot=FALSE, xlim=c(-2.3,2.3), template = "clock", 
    rotation = "clock", col="black", lwd = 3, lty=3, zero =c(rad(90)), main = " ",
    ylab = " ", xlab = " ", axes=FALSE, ticks=FALSE) # plot all
points(direction.rad[pds==1], rotation = "clock", zero =c(rad(90)), col = "black", pch=24 )
res1 <- density(direction.rad[cducsu==1], bw=5)
lines(res1, col="black", lwd = 2, lty=1, zero =c(rad(90)), rotation="clock")
points(direction.rad[cducsu==1], rotation = "clock", zero =c(rad(90)), col = "black", pch=22 )
res2 <- density(direction.rad[spd==1], bw=5, rotation="clock")
lines(res2, col="darkgrey", lwd = 2, lty=1, zero =c(rad(90)), rotation="clock")
points(direction.rad[spd==1], rotation = "clock", zero =c(rad(90)), col = "darkgrey", pch=25 )
res3 <- density(direction.rad[fdp==1], bw=5, rotation="clock")
lines(res3, col="gray", lwd = 3, lty=3, zero =c(rad(90)), rotation="clock")
points(direction.rad[fdp==1], rotation = "clock", zero =c(rad(90)), col = "gray", pch=21 )
res4 <- density(direction.rad[green==1], bw=5, rotation="clock")
lines(res4, col="darkgray", lwd = 3, lty=2, zero =c(rad(90)), rotation = "clock")
points(direction.rad[green==1], rotation = "clock", zero =c(rad(90)), 
    col = "darkgray", pch=23 )
legend("topright",c("CDU/CSU", "SPD", "FDP","Greens","PDS"),
    pch=c(22,25,21,23,24),
    col=c("black","darkgrey","gray","darkgray","black"), cex=1)
legend("topleft",c("CDU/CSU", "SPD", "FDP","Greens","PDS"),
    lty=c(1,1,3,2,3),lwd=c(2,2,3,3,3),
    col=c("black","darkgrey","gray","darkgray","black"), cex=1)
text(0,0.8,"0"); text(0,-0.8,expression(pi)); text(0.8,0,expression(pi/2)); text(-0.8,0,expression(3/2 *pi))
text(.5,-.1,"Dimension 1", cex=.5); text(0,.65,"Dimension 2", cex=.5); 
arrows(0, 0, 0, .5, angle= 0, code=1) 
arrows(0, 0, 0, -.5, angle= 0, code=1) 
arrows(0, 0, .5, 0, angle= 0, code=1) 
arrows(0, 0, -.5, 0, angle= 0, code=2) 

# COLOR
#### graph density german bundestag
set.seed(668)
res5 <- density(direction.rad[pds==1], bw=5)
plot(res5, points.plot=FALSE, xlim=c(-2.3,2.3), template = "clock", 
    rotation = "clock", col="darkred", zero =c(rad(90)), main = " ",
    ylab = " ", xlab = " ", axes=FALSE, ticks=FALSE) # plot all
points(direction.rad[pds==1], rotation = "clock", zero =c(rad(90)), col = "darkred", pch=24 )
res1 <- density(direction.rad[cducsu==1], bw=5)
lines(res1, col="black", zero =c(rad(90)), rotation="clock")
points(direction.rad[cducsu==1], rotation = "clock", zero =c(rad(90)), col = "black", pch=22 )
res2 <- density(direction.rad[spd==1], bw=5, rotation="clock")
lines(res2, col="red", zero =c(rad(90)), rotation="clock")
points(direction.rad[spd==1], rotation = "clock", zero =c(rad(90)), col = "red", pch=25 )
res3 <- density(direction.rad[fdp==1], bw=5, rotation="clock")
lines(res3, col="yellow3", zero =c(rad(90)), rotation="clock")
points(direction.rad[fdp==1], rotation = "clock", zero =c(rad(90)), col = "yellow3", pch=21 )
res4 <- density(direction.rad[green==1], bw=5, rotation="clock")
lines(res4, col="green", zero =c(rad(90)), rotation = "clock")
points(direction.rad[green==1], rotation = "clock", zero =c(rad(90)), 
    col = "green", pch=23 )
legend("topleft",c("CDU/CSU", "SPD", "FDP","Greens","PDS"),
    pch=c(22,25,21,23,24),
    col=c("black","red","yellow3","green","darkred"), cex=1)
text(0,0.8,"0"); text(0,-0.8,expression(pi)); text(0.8,0,expression(pi/2)); text(-0.8,0,expression(3/2 *pi))
text(.5,-.1,"Dimension 1", cex=.6); text(0,.65,"Dimension 2", cex=.6); 
arrows(0, 0, 0, .5, angle= 0, code=1) 
arrows(0, 0, 0, -.5, angle= 0, code=1) 
arrows(0, 0, .5, 0, angle= 0, code=1) 
arrows(0, 0, -.5, 0, angle= 0, code=2) 


source("C:/H/circular data/circ.lin.reg.R")

#### ANALYSIS
N <- length(year)
Y <- direction.rad
# APPROXIMATE CONSTANT ONLY-MODEL
out.null<-circ.lin.reg(matrix(1,N,1),Y, trace=T, tol = 1e-10, maxiter=100)


year6 <- year-min(year)
year6sq <- (year2^2)/100



unemp100 <- unemp/100
divorce100 <- divorce
outofwed100 <-outofwed/100
outofwed2.100 <-outofwed2/100


X <- cbind(year6,year6sq,unemp,outofwed2,reunification,spd,cducsu,green,pds)
out <- circ.lin.reg(X,Y,c(0), trace=T, tol = 1e-10, maxiter=1000)




# some numbers for Table 1 in the paper
out$log.lik # LL FULL MODEL
out.null$log.lik # LL NULL MODEL
(dev1 = 2 * (out$log.lik - out.null$log.lik)) # Fahrmeier/Tutz DEVIANCE
(dev2 = 2 * (-out$log.lik)) # RESIDUAL DEVIANCE
(dev3 = 2 * (0- out.null$log.lik)) # NULL DEVIANCE
(AIC <- 2*dim(X)[2] - 2*out$log.lik)  # AIC


out$mu / out$circ.se.mu 
out$kappa / out$se.kappa 

means <- mean(year6) * -0.169 + mean(year6sq) * 0.506 + mean(unemp) *  -0.084 + mean(outofwed2) *  -0.339
    + mean(reunification) * -0.501
#greens
(2 * atan(2.024 + means) - 2 * atan(means))
# PDS
(2 * atan(2.115 + means) - 2 * atan(means))
# SPD
(2 * atan(2.677 + means) - 2 * atan(means))
# CDU/CSU
(2 * atan(3.590 + means) - 2 * atan(means))
# FDP
(-0.8115243 + 2 * atan(0 + means))%%pi 

(2 * pi -0.8115243 + 2 * atan(0 + means))
# MH sampler for Fisher & Lee (1992)-style homoscedastic circular-linear regression model


library(circular)

ll.vm.atan <- function(b)  {
    N <- dim(X)[1]
    p <- dim(X)[2]+1
    coef  <- as.vector(b)
    kappa <- kappa
    mu <- coef[1]
    beta <- coef[2:p]
    ll1 <- - N * log(2 * pi * I.0(kappa))
    g.atan   <- mu + 2 * atan(X %*% beta) 
    ll2 <- kappa * sum(cos(Y  - g.atan ))
    ll <- ll1 + ll2 
    # sum priors  
    logpriors <- sum(dnorm(beta,0,10,log=TRUE)) + log(dvonmises(mu,out$mu,5))
    post <- ll + logpriors
    return(post)
}



### run "C:\H\circular data\germany\germany2.R here and save out <- circ.lin.reg


kappa <- out$kappa
kappa.var <- (out$se.kappa)^2
out$mu <- 2*pi - abs(out$mu)
out.par <- as.vector(c(out$mu,out$beta))
varcov.mat <- matrix(0,length(out$beta)+1,length(out$beta)+1)
varcov.mat[1,1] <- (out$circ.se.mu^2)/10  
varcov.mat[2:(length(out$beta)+1),2:(length(out$beta)+1)] <- out$cov.beta    
se <- sqrt(diag(abs(varcov.mat)))                     
z.val <- abs(out.par/se)
par.list <- cbind(out.par,as.vector(se),as.vector(z.val))
rownames(par.list) <- c("mu",colnames(cbind(X)))
colnames(par.list) <- c("coef", "se", "|z|")
print(round(par.list,3))
var.mat <- diag(1,dim(varcov.mat)[1]) * diag(varcov.mat)

mh.circ.lin=function(niter=1000,scale=.25,verbose=100)
{

# GAUSSIAN RANDOM WALK MH SAMPLER FOR THE HOMOSCEDASTIC CIRCULAR-LINEAR REGRESSION MODEL

    set.seed(666)
    library(mnormt)
    p=length(out.par)
    betas=matrix(0,niter,p)
    betas[1,]=as.vector(out.par)
    kappas=matrix(0,niter,1)
    kappas[1]=kappa
    likelihood=matrix(0,niter,1)
    likelihood[1]=ll.vm.atan(out.par) 
    Sigma2=as.matrix(var.mat)
    ver=seq(1,niter,by=verbose)
    
    for (i in 2:niter)
        {
        tildebeta=rmnorm(1,betas[i-1,],scale*Sigma2)
        tildebeta[1]= rvonmises(1, betas[i-1,1], 1/scale * 3)
        
        # GET TILDEKAPPA
            S <- sum(sin(Y - tildebeta[1] - 2 * atan(X %*% tildebeta[2:p])))/dim(X)[1]
            C <- sum(cos(Y - tildebeta[1]   - 2 * atan(X %*% tildebeta[2:p])))/dim(X)[1]
            R <- sqrt(S^2 + C^2)
            tildekappa=A1inv(R)
            kappa <- tildekappa
            llnew = ll.vm.atan(tildebeta)
            kappa <- kappas[i-1]
            llold = ll.vm.atan(betas[i-1,])
            llr = llnew - llold 

        draw=runif(1)
        if (draw<=exp(llr)) 
            {
            betas[i,]=tildebeta 
            kappas[i]=tildekappa
            likelihood[i]=llnew #+ llkappanew
            # get kappa
            S <- sum(sin(Y - betas[i,1] - 2 * atan(X %*% betas[i,2:p])))/dim(X)[1]
            C <- sum(cos(Y - betas[i,1] - 2 * atan(X %*% betas[i,2:p])))/dim(X)[1]
            R <- sqrt(S^2 + C^2)
            kappas[i]=A1inv(R)
            }
        else 
            {
            betas[i,]=betas[i-1,]
            kappas[i]=kappas[i-1]
            likelihood[i]=llold #+ llkappaold
            # get kappa
            S <- sum(sin(Y - betas[i,1] - 2 * atan(X %*% betas[i,2:p])))/dim(X)[1]
            C <- sum(cos(Y - betas[i,1] - 2 * atan(X %*% betas[i,2:p])))/dim(X)[1]
            R <- sqrt(S^2 + C^2)
            kappas[i]=A1inv(R)
            }
        if (is.na(match(i,ver))==FALSE) print(c(iter=i, betas[i,], kappas[i], like=likelihood[i]))
        }
 list(betas=betas, likelihood=likelihood, kappas=kappas)
}

library(mnormt)
post.samp=mh.circ.lin(niter=5000000, scale=.01, verbose=10000)
(ac.rate=accept(post.samp$betas)) # [1] 0.306
library(arm)
plot(as.mcmc(post.samp$betas))
plot(as.mcmc(post.samp$betas[,1]))

summary(as.mcmc(post.samp$betas))
summary(as.mcmc(post.samp$kappas))
plot(as.mcmc(post.samp$kappas))

# SAVE ALL DRAWS
#post.samp.full <- post.samp

plot(as.mcmc(post.samp$betas[,1:3]))
plot(as.mcmc(post.samp$betas[,4:6]))
plot(as.mcmc(post.samp$betas[,7:9]))

library(coda)
geweke.diag(post.samp$betas, frac1=0.1, frac2=0.5)
heidel.diag(post.samp$betas, eps=0.1, pvalue=0.05)


thin <- seq(1000000,5000000,4000)
post.samp$betas <- post.samp$betas[thin,]
post.samp$kappas<- post.samp$kappas[thin,]
post.samp$likelihood <- post.samp$likelihood[thin,]


# DEVIANCE INFORMATION CRITERION

# D.bar
(D.bar = -2*mean(post.samp$likelihood))
# D.theta
    ll1 <- - N * log(2 * pi * I.0(mean(post.samp$kappas)))
    g.atan   <- mean(post.samp$betas[,1]) + 2 * atan(X %*% apply(post.samp$betas[,2:(dim(X)[2]+1)],2,mean)) 
    ll2 <- mean(post.samp$kappas) * sum(cos(Y  - g.atan ))
(D.theta = -2*(ll1+ll2))
#P.D
(p.D = D.bar - D.theta)
#DIC
(DIC <- p.D + D.bar)


accept <- function(post.samp)
{
    t <- 0
    for(i in 2:dim(post.samp)[1])
    {
     if (post.samp[i,2]==post.samp[i-1,2])
    {
     t <- t + 1
     accept <- 1 - (t/dim(post.samp)[1])
    }
    }
   return(accept)
}
